
rm(list=ls())

library(tidyverse)
library(readxl)
library(here)
library(estimatr)
library(modelsummary)
library(here)
library(haven)
library(glmnet)
library(Matrix)
library(wesanderson)
library(ggpubr)
library(readxl)
library(fastDummies)
library(gridExtra)




########################## EXPERIMENT 1 - Informational experiment ################################## 

########

sl <- sl %>%
  mutate(treat = ifelse(is.na(FL_52_DO_MattB) | FL_52_DO_MattB != 1, 0, 1))
sl <- sl %>%
  mutate(treat = ifelse(is.na(treat), 0, treat))
table(sl$treat)

###Prepare the database ##### 

##### DEPENDENT VARIABLES: IMM1, IMM2, IMM3, IMM4 => INDEX

## IMM1: Ritiene che il numero di immigrati provenienti da paesi stranieri 
# che hanno un regolare permesso di trasferirsi in Italia dovrebbe essere...

# Named vector with numerical codes for the survey responses in Italian
immigration_responses <- c(
  'Significativamente ridotto' = 1,
  'leggermente ridotto' = 2,
  'lasciato invariato' = 3,
  'leggermente aumentato' = 4,
  'significativamente aumentato' = 5
)


sl <- sl %>%
  mutate(IMM1 = as.integer(factor(`immigration 1`, levels = names(immigration_responses), labels = immigration_responses)))



## IMM2: In generale, direbbe che gli immigrati clandestini offrono un contributo positivo alla societ� italiana 
# o rappresentano un fardello? 

# Mapping of string responses to integers
immigration2_responses <- c(
  '0 - Principalmente un fardello' = 0,
  '1' = 1,
  '2' = 2,
  '3' = 3,
  '4' = 4,
  '5 - Principalmente un contributo positivo' = 5
)


sl <- sl %>%
  mutate(IMM2 = as.integer(factor(`Immigration 2`, levels = names(immigration2_responses))))


# IMM3 In generale, direbbe che gli immigrati regolari offrono un contributo positivo 
# alla societ� italiana o rappresentano un fardello? 

immigration3_responses <- c(
  '0 - Principalmente un fardello' = 0,
  '1' = 1,
  '2' = 2,
  '3' = 3,
  '4' = 4,
  '5 - Principalmente un contributo positivo' = 5
)


sl <- sl %>%
  mutate(IMM3 = as.integer(factor(`Immigration 3`, levels = names(immigration3_responses))))

# IMM4_ E' favorevole o sfavorevole alla previsione di un percorso per ottenere la cittadinanza 
# per gli immigrati clandestini gi� presenti in Italia?

# Mapping of string responses to integers for 'Immigration 4'
immigration4_responses <- c(
  'Totalmente sfavorevole' = 1,
  'Leggermente sfavorevole' = 2,
  'Neutrale' = 3,
  'Leggermente favorevole' = 4,
  'Totalmente favorevole' = 5
)


sl <- sl %>%
  mutate(IMM4 = as.integer(factor(`Immigration 4`, levels = names(immigration4_responses))))



#### create an index (standardised and scaled)

sl <- sl %>%
  rowwise() %>%
  mutate(immigration_index = sum(c_across(c(IMM1, IMM2, IMM3, IMM4)), na.rm = TRUE)) %>%
  ungroup()
sl <- sl %>%
  mutate(immigration_index_standardized = scale(immigration_index, center = TRUE, scale = TRUE))



sl <- sl %>%
  rowwise() %>%
  mutate(immigration_index = if (any(is.na(c(IMM1, IMM2, IMM3, IMM4)))) NA else sum(c(IMM1, IMM2, IMM3, IMM4), na.rm = TRUE)) %>%
  ungroup() %>%
  mutate(immigration_index_standardized = ifelse(is.na(immigration_index), NA, scale(immigration_index, center = TRUE, scale = TRUE)))


sl <- sl %>%
  mutate(imm_vs_debt_coded = ifelse(imm_vs_debt == 'Ridurre il livello di debito pubblico', 1,
                                    ifelse(imm_vs_debt == 'Limitare l`immigrazione', 0, NA)))

sl <- sl %>%
  mutate(across(c(IMM1, IMM2, IMM3, IMM4), ~scale(.x, center = TRUE, scale = TRUE), .names = "standardized_{.col}"))


### drop respondents who didn't complete the survey #### 

sl <- sl %>% drop_na(immigration_index_standardized)


#### dep var (robustness test) measuring support for debt reduction & austerity

## Il governo dovrebbe ridurre il debito
sl <- sl %>%
  mutate(pubdebtA_coded = case_when(
    pubdebtA == 'Completamente in disaccordo' ~ 0,
    pubdebtA == 'In disaccordo' ~ 1,
    pubdebtA == 'N� d'accordo n� in disaccordo' ~ 2,
    pubdebtA == 'D'accordo' ~ 3,
    pubdebtA == 'Completamente d'accordo' ~ 4
  ))


sl <- sl %>%
  mutate(pubdebtB_coded = case_when(
    pubdebtB == 'Completamente in disaccordo' ~ 0,
    pubdebtB == 'In disaccordo' ~ 1,
    pubdebtB == 'N� d'accordo n� in disaccordo' ~ 2,
    pubdebtB == 'D'accordo' ~ 3,
    pubdebtB == 'Completamente d'accordo' ~ 4
  ))


sl <- sl %>%
  mutate(pubdebtC_coded = case_when(
    pubdebtC == 'Completamente in disaccordo' ~ 0,
    pubdebtC == 'In disaccordo' ~ 1,
    pubdebtC == 'N� d'accordo n� in disaccordo' ~ 2,
    pubdebtC == 'D'accordo' ~ 3,
    pubdebtC == 'Completamente d'accordo' ~ 4
  ))

## create an index

sl$pubdebtindex <- sl$pubdebtA_coded+sl$pubdebtB_coded+sl$pubdebtC_coded

sl <- sl %>%
  mutate(pubdebtindex_std = scale(pubdebtindex, center = TRUE, scale = TRUE))

sl <- sl %>%
  mutate(pubdebtA_std = scale(pubdebtA_coded, center = TRUE, scale = TRUE))
sl <- sl %>%
  mutate(pubdebtB_std = scale(pubdebtB_coded, center = TRUE, scale = TRUE))
sl <- sl %>%
  mutate(pubdebtC_std = scale(pubdebtC_coded, center = TRUE, scale = TRUE))


#######  Covariates ######

# age and partisanship 
sl$LR <- as.numeric(sl$left_right)
sl$LRsq <- sl$LR^2
sl$age <- as.numeric(sl$age)

# education & income 
table(sl$education)
table(sl$income)

sl <- sl %>%
  mutate(education_numeric = case_when(
    education == "Nessun titolo di studi" ~ 1,
    education == "Scuola elementare" ~ 2,
    education == "Scuola media" ~ 3,
    education == "Qualifica professionale (2 - 3 anni)" ~ 4,
    education == "Diploma di maturit� (5 anni)" ~ 5,
    education == "Laurea o post-laurea" ~ 6,
    TRUE ~ NA_real_  # Handle any unexpected values
  ))


sl <- sl %>%
  mutate(income_numeric = case_when(
    income == "Fino a 1000 euro" ~ 1,
    income == "Da 1.001 a 1.500 euro" ~ 2,
    income == "Da 1.501 a 2.000 euro" ~ 3,
    income == "Da 2.001 a 2.500 euro" ~ 4,
    income == "Da 2.501 a 3.000 euro" ~ 5,
    income == "Da 3.001 a 3.500 euro" ~ 6,
    income == "Da 3.501 a 4.000 euro" ~ 7,
    income == "Da 4.001 a 4.500 euro" ~ 8,
    income == "Da 4.501 a 5.000 euro" ~ 9,
    income == "Oltre 5.001 euro" ~ 10,
    TRUE ~ NA_real_  # Handle any unexpected values
  ))

# region and gender 

sl <- dummy_cols(sl, select_columns = c("region"), remove_first_dummy = TRUE, remove_selected_columns = TRUE)

sl <- sl %>%
  rename(
    region_Emilia_Romagna = `region_Emilia Romagna`,
    region_Friuli_Venezia_Giulia = `region_Friuli-Venezia Giulia`,
    region_Trentino_Alto_Adige = `region_Trentino-Alto Adige`,
    region_Valle_dAosta = `region_Valle d'Aosta`
  )

sl<-sl %>% mutate(woman = ifelse(gender=="Femmina",1,0))

sl_A <- sl %>% filter(treat == 0)
sl_B <- sl %>% filter(treat == 1)
mean(sl_A$IMM4, na.rm = TRUE)
mean(sl_B$IMM4, na.rm = TRUE)

### REGRESSION ANALYSIS ######


##### LASSO #####

##### lasso rule using cross validation. 
get_nonzero_coefficients <- function(rhs, dv) {
  # Perform cross-validated LASSO to select the best lambda
  cv_lasso <- cv.glmnet(rhs, dv, alpha = 1)
  
  # Get the lambda value that gives minimum cross-validated error (lambda.min)
  lambda_min <- cv_lasso$lambda.min
  
  # Alternatively, you could use the 1-SE rule lambda (lambda.1se) if you want a more conservative choice
  # lambda_min <- cv_lasso$lambda.1se
  
  # Extract the coefficients for the chosen lambda
  coefmat <- as.matrix(coef(cv_lasso, s = lambda_min))
  coefmat <- as.data.frame(coefmat)
  
  # Filter for non-zero coefficients
  nonzero <- coefmat %>% filter(s1 > 0)
  
  # Get the names of the non-zero coefficient variables
  row_names_vector <- as.vector(row.names(nonzero))
  row_names_vector <- row_names_vector[-1]  # Remove the intercept
  
  return(row_names_vector)
}

## select RHS varialbes; 
rhs<- sl %>%  select(LR,LRsq,age,education_numeric,woman,starts_with("region_")) %>%
  select(-contains("_DO_"))
rhs$`region_Da quale regione proviene?`<-NULL               

rhs <- rhs %>%
  mutate(across(where(is.numeric), ~ifelse(is.na(.), mean(., na.rm = TRUE), .)))
na_count <- sum(is.na(rhs))
print(na_count)

## turn rhs into matrix for the lasso package
rhs<-as.matrix(rhs)
# turn DV in to matrix for reasons.
dv <- as.matrix(sl$immigration_index_standardized)
dv1 <- as.matrix(sl$standardized_IMM1)
dv2 <- as.matrix(sl$standardized_IMM2)
dv3 <- as.matrix(sl$standardized_IMM3)
dv4 <- as.matrix(sl$standardized_IMM4)


row_names_vector <- get_nonzero_coefficients(rhs, dv)
print(row_names_vector)
row_names_vector1 <- get_nonzero_coefficients(rhs, dv1)
row_names_vector2 <- get_nonzero_coefficients(rhs, dv2)
row_names_vector3 <- get_nonzero_coefficients(rhs, dv3)
row_names_vector4 <- get_nonzero_coefficients(rhs, dv4)


# models
formula_index <- as.formula(paste("immigration_index_standardized ~ treat + ", paste(row_names_vector, collapse = "+")))
formula_IMM1 <- as.formula(paste("standardized_IMM1 ~ treat + ", paste(row_names_vector1, collapse = "+")))
formula_IMM2 <- as.formula(paste("standardized_IMM2 ~ treat + ", paste(row_names_vector2, collapse = "+")))
formula_IMM3 <- as.formula(paste("standardized_IMM3 ~ treat + ", paste(row_names_vector3, collapse = "+")))
formula_IMM4 <- as.formula(paste("standardized_IMM4 ~ treat + ", paste(row_names_vector4, collapse = "+")))


#plot
table1 <- list(
  "Index DV" = m1 <- lm_robust(formula_index, d=sl),
  "Support for increasing \n regular immigration" = m2 <- lm_robust(formula_IMM1, d=sl),
  "View of irregular immigrants" = m3 <- lm_robust(formula_IMM2, d=sl),
  "View of regular immigrants" = m4 <- lm_robust(formula_IMM3, d=sl),
  "Support for regularisation" = m5 <- lm_robust(formula_IMM4, d=sl)
  
)

modelsummary(table1, stars = TRUE)

modelsummary(table1, stars = TRUE, output = "maintable.docx" )

cmap<- ("Info. Treat" = "treat")

modelplot(table1, coef_map=cmap)


# Plot models and flip legend order, and remove y-axis label and text


finalgraph <- modelplot(table1, coef_map = cmap) +
  scale_color_manual(values = wes_palette("Darjeeling2")) +
  guides(color = guide_legend(reverse = TRUE)) + 
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey") +
  labs(y = NULL) +  # Remove y-axis label
  theme(axis.text.y = element_blank(),  # Remove y-axis text
        axis.ticks.y = element_blank()) +  # Remove y-axis ticks
  coord_cartesian(clip = 'off')  # Prevent clipping of elements

finalgraph

# Save the plot using ggsave, adjust size to remove extra space
ggsave("final_infortreat.png", finalgraph, width = 6, height = 4)


ggsave(
  filename = "final_infortreat.tiff",
  plot = finalgraph,
  device = "tiff",
  width = 6,
  height = 4,
  dpi = 300
)


############### CONDITIONAL EFFECTS ###################

sl <- sl %>%
  mutate(left = ifelse(LR <= 4, 0, 1))

###### conditional on LR 

formula_index_lr <- as.formula(paste("immigration_index_standardized ~ treat*left + ", paste(row_names_vector, collapse = "+")))
formula_IMM1_lr <- as.formula(paste("standardized_IMM1 ~ treat*left + ", paste(row_names_vector1, collapse = "+")))
formula_IMM2_lr <- as.formula(paste("standardized_IMM2 ~ treat*left + ", paste(row_names_vector2, collapse = "+")))
formula_IMM3_lr <- as.formula(paste("standardized_IMM3 ~ treat*left + ", paste(row_names_vector3, collapse = "+")))
formula_IMM4_lr <- as.formula(paste("standardized_IMM4 ~ treat*left + ", paste(row_names_vector4, collapse = "+")))


LRtable <- list(
  "Index DV" = m1_lr <- lm_robust(formula_index_lr, d=sl),
  "Support for increasing \n regular immigration" = m2_lr <- lm_robust(formula_IMM1_lr, d=sl),
  "View of irregular immigrants" = m3_lr <- lm_robust(formula_IMM2_lr, d=sl),
  "View of regular immigrants" = m4_lr <- lm_robust(formula_IMM3_lr, d=sl),
  "Support for regularisation" = m5_lr <- lm_robust(formula_IMM4_lr, d=sl)
  
)
modelsummary(LRtable, stars = TRUE)

modelsummary(LRtable, stars = TRUE, output = "LRtable.csv" )

## conditional on education 

sl <- sl %>% mutate(edu = education)

sl$edu[sl$edu == "Nessun titolo di studi" ] <- 1
sl$edu[sl$edu == "Scuola elementare" ] <- 1 
sl$edu[sl$edu == "Scuola media" ] <- 1
sl$edu[sl$edu == "Qualifica professionale (2 - 3 anni)" ] <- 2
sl$edu[sl$edu == "Diploma di maturit� (5 anni)" ] <- 2
sl$edu[sl$edu == "Laurea o post-laurea" ] <- 3

sl$edu <- as.factor(sl$edu)

formula_index_edu <- as.formula(paste("immigration_index_standardized ~ treat*edu + ", paste(row_names_vector, collapse = "+")))
formula_IMM1_edu <- as.formula(paste("standardized_IMM1 ~ treat*edu + ", paste(row_names_vector1, collapse = "+")))
formula_IMM2_edu <- as.formula(paste("standardized_IMM2 ~ treat*edu + ", paste(row_names_vector2, collapse = "+")))
formula_IMM3_edu <- as.formula(paste("standardized_IMM3 ~ treat*edu + ", paste(row_names_vector3, collapse = "+")))
formula_IMM4_edu <- as.formula(paste("standardized_IMM4 ~ treat*edu + ", paste(row_names_vector4, collapse = "+")))


EDUtable <- list(
  "Index DV" = m1_edu <- lm_robust(formula_index_edu, d=sl),
  "Support for increasing \n regular immigration" = m2_edu <- lm_robust(formula_IMM1_edu, d=sl),
  "View of irregular immigrants" = m3_edu <- lm_robust(formula_IMM2_edu, d=sl),
  "View of regular immigrants" = m4_edu <- lm_robust(formula_IMM3_edu, d=sl),
  "Support for regularisation" = m5_edu <- lm_robust(formula_IMM4_edu, d=sl)
  
)
modelsummary(EDUtable, stars = TRUE)

modelsummary(EDUtable, stars = TRUE, output = "EDU_table.docx" )



######### ROBUSTNESS checks #### 

### austerity outcome: 

sl <- sl %>% drop_na(pubdebtindex_std)

## select RHS varialbes; 
rhs<- sl %>%  select(LR,LRsq,age,education_numeric,woman,starts_with("region_")) %>%
  select(-contains("_DO_"))
rhs$`region_Da quale regione proviene?`<-NULL               

rhs <- rhs %>%
  mutate(across(where(is.numeric), ~ifelse(is.na(.), mean(., na.rm = TRUE), .)))
na_count <- sum(is.na(rhs))
print(na_count)

## turn rhs into matrix for the lasso package
rhs<-as.matrix(rhs)
# turn DV in to matrix for reasons.
dv <- as.matrix(sl$pubdebtindex_std)
dv1 <- as.matrix(sl$pubdebtA_std)
dv2 <- as.matrix(sl$pubdebtB_std)
dv3 <- as.matrix(sl$pubdebtC_std)


row_names_vector <- get_nonzero_coefficients(rhs, dv)
print(row_names_vector)
row_names_vector1 <- get_nonzero_coefficients(rhs, dv1)
row_names_vector2 <- get_nonzero_coefficients(rhs, dv2)
row_names_vector3 <- get_nonzero_coefficients(rhs, dv3)


formula_index_aus <- if (length(row_names_vector) > 0) {
  as.formula(paste("pubdebtindex_std ~ treat +", paste(row_names_vector, collapse = "+")))
} else {
  as.formula("pubdebtindex_std ~ treat")
}

formula_ausA <- if (length(row_names_vector1) > 0) {
  as.formula(paste("pubdebtA_std ~ treat +", paste(row_names_vector1, collapse = "+")))
} else {
  as.formula("pubdebtA_std ~ treat")
}

formula_ausB <- if (length(row_names_vector2) > 0) {
  as.formula(paste("pubdebtB_std ~ treat +", paste(row_names_vector2, collapse = "+")))
} else {
  as.formula("pubdebtB_std ~ treat")
}

formula_ausC <- if (length(row_names_vector3) > 0) {
  as.formula(paste("pubdebtC_std ~ treat +", paste(row_names_vector3, collapse = "+")))
} else {
  as.formula("pubdebtC_std ~ treat")
}

##tablr 
table1 <- list(
  "Austerity Index" = m1_aus <- lm_robust(formula_index_aus, d=sl),
  "No trade-off" = m2_aus <- lm_robust(formula_ausA, d=sl),
  "Spending Cuts" = m3_aus <- lm_robust(formula_ausB, d=sl),
  "Taxes" = m4_aus <- lm_robust(formula_ausC, d=sl)
)

modelsummary(table1)
modelsummary(table1, stars = TRUE, output = "austeritytable.csv" )

##### Figure 4 main text ##### 

# Plot austerity models

finalgraph_aus <- modelplot(table1, coef_map = cmap) +
  scale_color_manual(values = wes_palette("Darjeeling2")) +
  guides(color = guide_legend(reverse = TRUE)) + 
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey") +
  labs(y = NULL) +  # Remove y-axis label
  theme(axis.text.y = element_blank(),  # Remove y-axis text
        axis.ticks.y = element_blank())  # Remove y-axis ticks


finalgraph_aus

# Save the plot using ggsave, adjust size to remove extra space
ggsave(filename = "finalgraph_aus.tiff",
       plot = finalgraph_aus,
       device = "tiff",
       width = 6,
       height = 4,
       dpi = 300)

############### attention check 

sl <- sl %>% 
  mutate(passcheck = if_else(att_check == "Altri quotidiani", 1, 0))
table(sl$passcheck)

#subset
sl_att <- sl %>% filter(passcheck == 1)


## select RHS varialbes
rhs<- sl_att %>%  select(LR,LRsq,age,education_numeric,woman,starts_with("region_")) %>%
  select(-contains("_DO_"))
rhs$`region_Da quale regione proviene?`<-NULL               

rhs <- rhs %>%
  mutate(across(where(is.numeric), ~ifelse(is.na(.), mean(., na.rm = TRUE), .)))
na_count <- sum(is.na(rhs))
print(na_count)

## turn rhs into matrix for the lasso package
rhs<-as.matrix(rhs)

# turn DV in to matrix
dv_att <- as.matrix(sl_att$immigration_index_standardized)
dv1_att <- as.matrix(sl_att$standardized_IMM1)
dv2_att <- as.matrix(sl_att$standardized_IMM2)
dv3_att <- as.matrix(sl_att$standardized_IMM3)
dv4_att <- as.matrix(sl_att$standardized_IMM4)


row_names_vector <- get_nonzero_coefficients(rhs, dv_att)
print(row_names_vector)
row_names_vector1 <- get_nonzero_coefficients(rhs, dv1_att)
row_names_vector2 <- get_nonzero_coefficients(rhs, dv2_att)
row_names_vector3 <- get_nonzero_coefficients(rhs, dv3_att)
row_names_vector4 <- get_nonzero_coefficients(rhs, dv4_att)


# models
formula_index <- as.formula(paste("immigration_index_standardized ~ treat + ", paste(row_names_vector, collapse = "+")))
formula_IMM1 <- as.formula(paste("standardized_IMM1 ~ treat + ", paste(row_names_vector1, collapse = "+")))
formula_IMM2 <- as.formula(paste("standardized_IMM2 ~ treat + ", paste(row_names_vector2, collapse = "+")))
formula_IMM3 <- as.formula(paste("standardized_IMM3 ~ treat + ", paste(row_names_vector3, collapse = "+")))
formula_IMM4 <- as.formula(paste("standardized_IMM4 ~ treat + ", paste(row_names_vector4, collapse = "+")))


#plot
table_att <- list(
  "Index DV" = m1_att <- lm_robust(formula_index, d=sl_att),
  "Support for increasing \n regular immigration" = m2_att <- lm_robust(formula_IMM1, d=sl_att),
  "View of irregular immigrants" = m3_att <- lm_robust(formula_IMM2, d=sl_att),
  "View of regular immigrants" = m4_att <- lm_robust(formula_IMM3, d=sl_att),
  "Support for regularisation" = m5_att <- lm_robust(formula_IMM4, d=sl_att)
  
)

modelsummary(table_att, stars = TRUE)

modelsummary(table_att, stars = TRUE, output = "maintable_att.docx" )



########### END ############################## 





###################### EXPERIMENT 2: Conjoint experiment ###################################�

library(gridExtra)
library(grDevices)
library(ggplot2)
library(tidyverse)
library(cjoint)
library(jtools)
library(interplot)
library(cregg)
library(gtable)
library(fastDummies)
library(modelsummary)
library(lmtest)
library(sandwich)
library(car)





setwd() ######## IMPORTANT: before starting, set the working director where the trilemma.csv file is stored 
#### the covariates have to be in the same order as they are in the csv file 

# check always for the automatic "respondent" column 

trilemma_db <- read.qualtrics("trilemma.csv", 
                              # the DV for our conjoint tasks (Candidate 1 and Candidate 2)
                              responses=c("task1_0", "task2_0","task3_0", "task4_0", "task5_0", "task6_0", 
                                          "task7_0", "task8_0"), 
                              covariates=c( "age", "gender", "education", "income", "left_right", "region",
                                            "t1_p1", "t1_p2", 
                                            "t2_p1", "t2_p2", "t3_p1", "t3_p2", "t4_p1", "t4_p2", "t5_p1", "t5_p2",
                                            "t6_p1", "t6_p2", "t7_p1", "t7_p2", "t8_p1", "t8_p2", "att_check", "party"),
                              respondentID = "ResponseId", letter = "F",  new.format= TRUE, ranks = NULL)


# remove respondents who didn't complete the survey 
trilemma_db <- trilemma_db %>% drop_na(selected) 

##################### 



### Measure the AMCE of pro-immigration policies vs broad-based austerity 

# create a new variable for the debt reduction strategy attribute, called "Debt_reduction" 
## and collapse pro-immigration policies on the one hand and broad-based austerity measures on the other 

trilemma_db <- trilemma_db %>% mutate(Debt_reduction = ifelse(Riduzione.del.debito == "Agevolare l'immigrazione di extracomunitari altamente qualificati"
                                                              |Riduzione.del.debito == "Agevolare l'immigrazione di Europei altamente qualificati"
                                                              |Riduzione.del.debito == "Agevolare l'immigrazione di Europei anche non qualificati"
                                                              |Riduzione.del.debito == "Agevolare l'immigrazione di extracomunitari anche non qualificati",
                                                              "Pro-Immigration Policies", 
                                                              ifelse(Riduzione.del.debito == "Aumento tasse sui redditi maggiori di 45.000 euro lordi annui", "High-income tax",
                                                                     ifelse(Riduzione.del.debito == "Nessuna misura per ridurre il debito", "No debt reduction",
                                                                            "Broad-based austerity"
                                                                            
                                                                     ))))

### compute the AMCE (Figure 5) main text ######### 

trilemma_db$Debt_reduction <- as.factor(trilemma_db$Debt_reduction)

# Change the baseline for fiscal policies in "broad-based austerity"

baselines <- list()
baselines$Debt_reduction <- "Broad-based austerity"


eng_coll <- cjoint::amce(selected ~ Debt_reduction,
                         data= trilemma_db,
                         cluster=TRUE, 
                         respondent.id= "Response.ID", baselines = baselines)


plot(eng_coll, attribute.names = "Debt reduction strategies",
     xlab = "Estimated AMCE",
     plot.theme = theme(text = element_text(size = 14,
                                            face = c("plain", "plain", "plain","italic","bold")),
                        legend.position = 'none',
                        panel.grid.major= element_blank(),
                        panel.grid.minor= element_blank(),
                        panel.background = element_blank(),axis.line = element_line(colour = "black")
     ))


summary(eng_coll) 

# Open a TIFF graphics device
tiff("debt_reduction_strategies.tiff", width = 6, height = 4, units = "in", res = 300)

# Draw the plot
plot(eng_coll,
     attribute.names = "Debt reduction strategies",
     xlab = "Estimated AMCE",
     plot.theme = theme(
       text = element_text(size = 14,
                           face = c("plain", "plain", "plain", "italic", "bold")),
       legend.position = 'none',
       panel.grid.major = element_blank(),
       panel.grid.minor = element_blank(),
       panel.background = element_blank(),
       axis.line = element_line(colour = "black")
     )
)

# Close the device
dev.off()







#### Further analysis with mm (supplementary material, appendix C) ######### 



mm_details <- mm(trilemma_db, selected ~ Riduzione.del.debito,
                 cluster=TRUE, 
                 respondent.id= "Response.ID")

custom_labels <- c( "Pro-immigration policies - EU highly skilled",
                    "Pro-immigration policies - EU unskilled",
                    "Pro-immigration policies - non-EU highly skilled",
                    "Pro-immigration policies - non-EU unskilled",
                    "Tax increases - VAT",
                    "Tax increases - Income above 45k",
                    "No debt reduction strategy",
                    "Spending cuts - education & childcare",
                    "Spending cuts - pensions & unemployment benefits")


mm_plot <- ggplot(mm_details, aes(y = level, x = estimate)) +
  geom_point() +
  geom_errorbar(aes(xmin = lower, xmax = upper), width = 0, linetype = "solid") +
  geom_vline(xintercept = 0.5, linetype = "dashed", color = "black") +
  labs(title = "",
       x = "Estimated Marginal Means",
       y = "") +
  scale_y_discrete(labels = custom_labels) +
  theme_minimal()

summary(mm_details)





################### MMs with all of the attributes (supplementary material C) #######

trilemma_db <- trilemma_db %>% mutate(Defense = Difesa)
trilemma_db <- trilemma_db %>% mutate(Civil_rights = Diritti.e.tutela.dei.cittadini)
trilemma_db <- trilemma_db %>% mutate(Energy_independence  = Indipendenza.energetica)

#rename attributes and levels 

# Defense

trilemma_db$Defense <- as.character(trilemma_db$Defense)

trilemma_db$Defense[grepl("Aumentare", trilemma_db$Defense)] <- "EU military capacity and single European army"
trilemma_db$Defense[grepl("Promozione di", trilemma_db$Defense)] <- "International initiatives for multilateral disarmament"
trilemma_db$Defense[grepl("Rafforzamento", trilemma_db$Defense)] <- "European diplomacy"
trilemma_db$Defense[grepl("Sanzioni", trilemma_db$Defense)] <- "Sanctions for selling weapons to countries at war"
trilemma_db$Defense[grepl("Mantere lo", trilemma_db$Defense)] <- "Status quo and NATO"

trilemma_db$Defense <- as.factor(trilemma_db$Defense)

#Energy

trilemma_db$Energy_independence <- as.character(trilemma_db$Energy_independence)

trilemma_db$Energy_independence[grepl("Efficientamento energetico", trilemma_db$Energy_independence)] <- "Energy efficiency and electrified public transport "
trilemma_db$Energy_independence[grepl("Diversificazione fonti", trilemma_db$Energy_independence)] <- "Diversified sources including nuclear energy"
trilemma_db$Energy_independence[grepl("Unione energetica", trilemma_db$Energy_independence)] <- "Energy Union"
trilemma_db$Energy_independence[grepl("Supporto alle", trilemma_db$Energy_independence)] <- "Energy communities for renewable energy"

trilemma_db$Energy_independence <- as.factor(trilemma_db$Energy_independence)


# Rights 
trilemma_db$Civil_rights <- as.character(trilemma_db$Civil_rights)

trilemma_db$Civil_rights[grepl("Rafforzare la", trilemma_db$Civil_rights)] <- "European legislation on gender violence"
trilemma_db$Civil_rights[grepl("Garantire la", trilemma_db$Civil_rights)] <- "Equal rights for same-sex couples"
trilemma_db$Civil_rights[grepl("Introduzione del diritto", trilemma_db$Civil_rights)] <- "Right to assisted death"
trilemma_db$Civil_rights[grepl("Quadro europeo per garantire", trilemma_db$Civil_rights)] <- "European framework for dignified detention"
trilemma_db$Civil_rights[grepl("Misure a tutela", trilemma_db$Civil_rights)] <- "Protection of independence of journalism"
trilemma_db$Civil_rights[grepl("Inserire il", trilemma_db$Civil_rights)] <- "Right to abortion in the Constitution"

trilemma_db$Civil_rights <- as.factor(trilemma_db$Civil_rights)

# mm 
full_mm <- mm(selected ~ Debt_reduction + Defense +  Civil_rights + Energy_independence,
                          data=trilemma_db,
                          cluster=TRUE, 
                          respondent.id= "Response.ID", baselines = baselines)


plot(full_mm)  +
  geom_vline(xintercept = 0.5, linetype = "dotted", color = "black")


####### EXPLORATORY analysis #########################################

############ SKILLS AND COUNTRY OF ORIGIN ###########

trilemma_db <- trilemma_db %>% mutate(Debt_reduction_2 = Riduzione.del.debito)

##### preliminary: add the new levels to the factor variable Debt_reduction_2

#imm
trilemma_db$Debt_reduction_2 <- factor(trilemma_db$Debt_reduction_2, levels = c(levels(trilemma_db$Debt_reduction_2), "Pro-immigration policies - EU highly skilled"))
trilemma_db$Debt_reduction_2 <- factor(trilemma_db$Debt_reduction_2, levels = c(levels(trilemma_db$Debt_reduction_2), "Pro-immigration policies - non-EU highly skilled"))
trilemma_db$Debt_reduction_2 <- factor(trilemma_db$Debt_reduction_2, levels = c(levels(trilemma_db$Debt_reduction_2), "Pro-immigration policies - EU unskilled"))
trilemma_db$Debt_reduction_2 <- factor(trilemma_db$Debt_reduction_2, levels = c(levels(trilemma_db$Debt_reduction_2), "Pro-immigration policies - non-EU unskilled"))
#aus 

trilemma_db$Debt_reduction_2 <- factor(trilemma_db$Debt_reduction_2, levels = c(levels(trilemma_db$Debt_reduction_2), "Broad-based austerity"))

trilemma_db$Debt_reduction_2 <- factor(trilemma_db$Debt_reduction_2, levels = c(levels(trilemma_db$Debt_reduction_2), "No debt reduction strategy"))

trilemma_db$Debt_reduction_2 <- factor(trilemma_db$Debt_reduction_2, levels = c(levels(trilemma_db$Debt_reduction_2), "Increase taxes - Income above 45k"))


#####
#rename the immigration policies (do not group them together this time)
trilemma_db$Debt_reduction_2[trilemma_db$Debt_reduction_2 == "Agevolare l'immigrazione di Europei altamente qualificati" ] <- "Pro-immigration policies - EU highly skilled"   
trilemma_db$Debt_reduction_2[trilemma_db$Debt_reduction_2 == "Agevolare l'immigrazione di extracomunitari altamente qualificati" ] <- "Pro-immigration policies - non-EU highly skilled"   
trilemma_db$Debt_reduction_2[trilemma_db$Debt_reduction_2 == "Agevolare l'immigrazione di Europei anche non qualificati" ] <- "Pro-immigration policies - EU unskilled"   
trilemma_db$Debt_reduction_2[trilemma_db$Debt_reduction_2 == "Agevolare l'immigrazione di extracomunitari anche non qualificati" ] <- "Pro-immigration policies - non-EU unskilled"   

# rename the aus policies (group them)
trilemma_db$Debt_reduction_2[grepl("Taglio spesa pubblica in pensioni", trilemma_db$Debt_reduction_2)] <- "Broad-based austerity"
trilemma_db$Debt_reduction_2[grepl("Taglio della spesa pubblica in istruzione", trilemma_db$Debt_reduction_2)] <- "Broad-based austerity"
trilemma_db$Debt_reduction_2[trilemma_db$Debt_reduction_2 == "Aumento dell'IVA" ] <- "Broad-based austerity"   

trilemma_db$Debt_reduction_2[trilemma_db$Debt_reduction_2 == "Nessuna misura per ridurre il debito" ] <- "No debt reduction strategy"   
trilemma_db$Debt_reduction_2[trilemma_db$Debt_reduction_2 == "Aumento tasse sui redditi maggiori di 45.000 euro lordi annui" ] <- "Increase taxes - Income above 45k" 


### AMCE Figure 6 Main text ###

baselines_3 <- list()
baselines_3$Debt_reduction_2 <- "Broad-based austerity"

m1 <- cjoint::amce(selected ~ Debt_reduction_2,
                   data= trilemma_db,
                   cluster=TRUE, 
                   respondent.id= "Response.ID", baselines = baselines_3)


plot(m1, attribute.names = "Debt reduction strategies",
     xlab = "Estimated AMCE",
     plot.theme = theme(text = element_text(size = 14,
                                            face = c("plain","plain","plain","plain", "plain", "plain","italic","bold")),
                        legend.position = 'none',
                        panel.grid.major= element_blank(),
                        panel.grid.minor= element_blank(),
                        panel.background = element_blank(),axis.line = element_line(colour = "black")
     ))


summary(m1) 


### test whether the coefficients are different 


# lm for hypotheses testing

trilemma_db$Debt_reduction_2 <- relevel(trilemma_db$Debt_reduction_2, ref = "Broad-based austerity")

m1 <- lm(selected ~ Debt_reduction_2, data = trilemma_db)
summary(m1)


# hpotheses
# Ha- B1_imm_EU_hs > B2_imm_NonEU_hs 
# Hb- B3_imm_EU_ls > B4_imm_NonEU_ls ;
# Hc- B1_imm_EU_hs > B3_imm_EU_ls ;
# Hd - B2_imm_NonEU_hs > B4_imm_NonEU_ls 


B1_imm_EU_hs <- 0.174482
B2_imm_NonEU_hs <- 0.173742
B3_imm_EU_ls <- 0.124143
B4_imm_NonEU_ls <- 0.081244

# calculate the difference (id it's positive) and next, whether the difference is significant
B1_imm_EU_hs - B2_imm_NonEU_hs # positive but close to zero
B3_imm_EU_ls - B4_imm_NonEU_ls # positive
B1_imm_EU_hs - B3_imm_EU_ls #positive 
B2_imm_NonEU_hs - B4_imm_NonEU_ls #positive

#


# Test hypothesis a
linearHypothesis(m1, "Debt_reduction_2Pro-immigration policies - EU HIGLY qualified = Debt_reduction_2Pro-immigration policies - non-EU HIGHLY qualified")

"No significant difference, cannot reject the null, i.e., reject Ha"

# Test hypothesis b
linearHypothesis(m1, "Debt_reduction_2Pro-immigration policies - EU NOT qualified = Debt_reduction_2Pro-immigration policies - non-EU NOT qualified")

"The p-value is less than the conventional significance level of 0.05, 
=> there is a statistically significant difference as hypothesised by Hb; origin matters only for low skilled"

# Test hypothesis c
linearHypothesis(m1, "Debt_reduction_2Pro-immigration policies - EU HIGLY qualified  = Debt_reduction_2Pro-immigration policies - EU NOT qualified")

"Result: F(1,18729) = 11.036, ????=  0.0008954
=> There is a statistically significant difference between the effects of 
Pro-immigration policies - EU HIGHLY qualified (B1_imm_EU_hs)
and Pro-immigration policies - EU NOT qualified (B2_imm_EU_ls),
on the probability of selecting a profile, with the former having a stronger impact. Skills matter"

# Test hypothesis d
linearHypothesis(m1, "Debt_reduction_2Pro-immigration policies - non-EU HIGHLY qualified = Debt_reduction_2Pro-immigration policies - non-EU NOT qualified")

"Very strong statistical evidence to reject the null hypothesis. Skills matter"


############### DIFFERENT AUSTERITY MEASURES #############

#create a new variable
trilemma_db <- trilemma_db %>% mutate(Debt_reduction_3 = Riduzione.del.debito)

trilemma_db$Debt_reduction_3 <- factor(trilemma_db$Debt_reduction_3, levels = c(levels(trilemma_db$Debt_reduction_3), "Tax increases - Income above 45k"))
trilemma_db$Debt_reduction_3 <- factor(trilemma_db$Debt_reduction_3, levels = c(levels(trilemma_db$Debt_reduction_3), "Tax Increases - VAT"))
trilemma_db$Debt_reduction_3 <- factor(trilemma_db$Debt_reduction_3, levels = c(levels(trilemma_db$Debt_reduction_3), "Spending cuts - education & childcare"))
trilemma_db$Debt_reduction_3 <- factor(trilemma_db$Debt_reduction_3, levels = c(levels(trilemma_db$Debt_reduction_3), "Spending cuts - pensions & unemployment benefits"))

trilemma_db$Debt_reduction_3 <- factor(trilemma_db$Debt_reduction_3, levels = c(levels(trilemma_db$Debt_reduction_3), "Pro-immigration policies"))
trilemma_db$Debt_reduction_3 <- factor(trilemma_db$Debt_reduction_3, levels = c(levels(trilemma_db$Debt_reduction_3), "No debt reduction strategy"))



# rename the aus policies (DO NOT group them this time)
trilemma_db$Debt_reduction_3[grepl("Taglio spesa pubblica in pensioni", trilemma_db$Debt_reduction_3)] <- "Spending cuts - pensions & unemployment benefits"
trilemma_db$Debt_reduction_3[grepl("Taglio della spesa pubblica in istruzione", trilemma_db$Debt_reduction_3)] <- "Spending cuts - education & childcare"
trilemma_db$Debt_reduction_3[trilemma_db$Debt_reduction_3 == "Aumento dell'IVA" ] <- "Tax Increases - VAT"   
trilemma_db$Debt_reduction_3[grepl("Aumento tasse sui", trilemma_db$Debt_reduction_3)] <- "Tax increases - Income above 45k"

#
trilemma_db$Debt_reduction_3[trilemma_db$Debt_reduction_3 == "Nessuna misura per ridurre il debito" ] <- "No debt reduction strategy"   

#
#rename the immigration policies (group them together )
trilemma_db$Debt_reduction_3[trilemma_db$Debt_reduction_3 == "Agevolare l'immigrazione di Europei altamente qualificati" ] <- "Pro-immigration policies"   
trilemma_db$Debt_reduction_3[trilemma_db$Debt_reduction_3 == "Agevolare l'immigrazione di extracomunitari altamente qualificati" ] <- "Pro-immigration policies"   
trilemma_db$Debt_reduction_3[trilemma_db$Debt_reduction_3 == "Agevolare l'immigrazione di Europei anche non qualificati" ] <- "Pro-immigration policies"   
trilemma_db$Debt_reduction_3[trilemma_db$Debt_reduction_3 == "Agevolare l'immigrazione di extracomunitari anche non qualificati" ] <- "Pro-immigration policies"   


baselines_4 <- list()
baselines_4$Debt_reduction_3 <- "Pro-immigration policies"

### FIGURE 7 main text ###

m2 <- cjoint::amce(selected ~ Debt_reduction_3,
                   data= trilemma_db,
                   cluster=TRUE, 
                   respondent.id= "Response.ID", baselines = baselines_4)


plot(m2, attribute.names = "Debt reduction strategies",
     xlab = "Estimated AMCE",
     plot.theme = theme(text = element_text(size = 14,
                                            face = c("plain","plain","plain","plain", "plain", "italic","bold")),
                        legend.position = 'none',
                        panel.grid.major= element_blank(),
                        panel.grid.minor= element_blank(),
                        panel.background = element_blank(),axis.line = element_line(colour = "black")
     ))


summary(m2) 


### test whether the coefficients are different 


#lm 
trilemma_db$Debt_reduction_3 <- relevel(trilemma_db$Debt_reduction_3, ref = "No debt reduction strategy")


m2 <- lm(selected ~ Debt_reduction_3, data = trilemma_db)
summary(m2)



#tax-based aus
# Hypothesis: Debt_reduction_3Tax increases - Income above 45k = Debt_reduction_3Tax Increases - VAT
hypothesis <- "Debt_reduction_3Tax increases - Income above 45k = Debt_reduction_3Tax Increases - VAT"

# Perform the linear hypothesis test
test_result <- linearHypothesis(m2, hypothesis)


# Spending-based aus
hypothesis <- "Debt_reduction_3Spending cuts - pensions & unemployment benefits = Debt_reduction_3Spending cuts - education & childcare"

# Linear hypothesis test
test_result <- linearHypothesis(m2, hypothesis)

# end of exploratory analysis 

#### CONDITIONAL effects #######

# partisanship (Figure 8) main text

baselines <- list()
baselines$Debt_reduction <- "Broad-based austerity"

trilemma_db <- trilemma_db %>%
  filter(!is.na(left_right)) %>%  # Exclude NA values in 'left_right'
  mutate(Partisanship = ifelse(left_right <= 4, 0, 1)) %>%
  mutate(Partisanship = as.factor(Partisanship))

summary(eng_coll) 

x_limits <- c(0.35, 0.7) # Define the x-axis limits

# Calculate marginal means and plot
fig1 <- plot(cregg::cj(trilemma_db, selected ~ Debt_reduction, 
               id = ~Response.ID,
               by = ~Partisanship,
               estimate = "mm"), 
     group = "Partisanship") +
  scale_color_manual(values = c("0" = "red", "1" = "blue"), 
                     labels = c("0" = "Left", "1" = "Right"),
                     na.translate = FALSE) +
  geom_vline(xintercept = 0.5, linetype = "dotted", color = "black") +  
  coord_cartesian(xlim = x_limits) +
  labs(x = NULL) + 
  guides(color = guide_legend(reverse = TRUE)) + # Remove x-axis title
  theme_bw() +
  theme(axis.text.y = element_text(size = 12)) 



#### education (Figure 8) main text

trilemma_db <- trilemma_db %>% mutate(Education = education)

trilemma_db$Education[trilemma_db$Education == 1 ] <- 1
trilemma_db$Education[trilemma_db$Education == 2 ] <- 1
trilemma_db$Education[trilemma_db$Education == 3 ] <- 1
trilemma_db$Education[trilemma_db$Education == 4 ] <- 2
trilemma_db$Education[trilemma_db$Education == 5 ] <- 2
trilemma_db$Education[trilemma_db$Education == 6 ] <- 3

baselines <- list()
baselines$Debt_reduction <- "Broad-based austerity"

trilemma_db$Education <- as.factor(trilemma_db$Education)


fig2 <- plot(cregg::cj(trilemma_db, selected ~ Debt_reduction, 
               id = ~Response.ID,
               by = ~Education,
               estimate = "mm"), 
     group = "Education") +
  scale_color_manual(values = c("1" = "#d95f02", "2" = "#1b9e77", "3" = "#7570b3"), 
                     labels = c("1" = "Education 1", "2" = "Education 2", "3" = "Education 3"),
                     na.translate = FALSE) +  # Remove NA from legend)+
  geom_vline(xintercept = 0.5, linetype = "dotted", color = "black") +  
  coord_cartesian(xlim = x_limits) +
  labs(x = NULL) + 
guides(color = guide_legend(reverse = TRUE)) + # Remove x-axis title
   theme_bw()  +
  theme(axis.text.y = element_text(size = 12)) 

grid.arrange(fig1, fig2, ncol = 2, bottom = "Marginal Means")
  

# respondents' region of residence (Appendix F)

trilemma_db <- trilemma_db %>%
  mutate(macroregion = case_when(
    region %in% c(1, 2, 3, 4, 11, 12) ~ "South",
    region %in% c(5, 6, 17, 20) ~ "North-Est",
    region %in% c(7, 10, 16, 18) ~ "Centre",
    region %in% c(14, 15) ~ "Islands",
    region %in% c(8, 9, 12,19) ~ "North-West",
    TRUE ~ "Unknown"
  ))

trilemma_db <- trilemma_db %>% filter(macroregion != "Unknown")

trilemma_db$macroregion <- as.factor(trilemma_db$macroregion)

m_region <- lm(selected ~ Debt_reduction*macroregion, data = trilemma_db )

m_region <- cregg::cj(trilemma_db, selected ~ Debt_reduction, 
                id = ~ Response.ID,
                by = ~ macroregion,
                estimate = "mm")

summary(m_region)

trilemma_db$macroregion <- as.factor(trilemma_db$macroregion)

region_fig <- plot(cregg::cj(trilemma_db, selected ~ Debt_reduction, 
                             id = ~Response.ID,
                             by = ~macroregion,
                             estimate = "mm"), 
                   group = "macroregion") +
  geom_vline(xintercept = 0.5, linetype = "dotted", color = "black") 



region_fig




### CENSUS, Appendix D ####### 

## Exp 1

# age

sl$age <- as.numeric(sl$age)

age_groups <- sl %>%
  mutate(age_group = case_when(
    age >= 18 & age <= 29 ~ "18-29",
    age >= 30 & age <= 39 ~ "30-39",
    age >= 40 & age <= 49 ~ "40-49",
    age >= 50 & age <= 59 ~ "50-59",
    age >= 60 ~ "60+",
    TRUE ~ NA_character_)) %>%
  group_by(age_group) %>%
  summarise(
    count = n(),
    percentage = (n() / nrow(sl)) * 100
  ) %>%
  filter(!is.na(age_group))  # Exclude any NA groups

print(age_groups)  

# regions

sl <- sl %>%
  mutate(macroregion = case_when(
    residence %in% c(1, 2, 3, 4, 11, 14) ~ "Sud",
    residence %in% c(5, 6, 18, 21) ~ "Nord-Est",
    residence %in% c(7, 10, 17, 19) ~ "Centro",
    residence %in% c(15, 16) ~ "Isole",
    residence %in% c(8, 9, 12, 20) ~ "Nord-Ovest",
    TRUE ~ "Unknown"
  ))

macroregion_summary <- sl %>%
  group_by(macroregion) %>%
  summarise(
    count = n(),
    percentage = (n() / nrow(sl)) * 100
  ) %>%
  filter(macroregion != "Unknown")  # Exclude any unknown macroregions

print(macroregion_summary)  

## exp 2

#age

age_groups <- trilemma_db %>%
  mutate(age_group = case_when(
    age >= 18 & age <= 29 ~ "18-29",
    age >= 30 & age <= 39 ~ "30-39",
    age >= 40 & age <= 49 ~ "40-49",
    age >= 50 & age <= 59 ~ "50-59",
    age >= 60 ~ "60+",
    TRUE ~ NA_character_)) %>%
  group_by(age_group) %>%
  summarise(
    count = n(),
    percentage = (n() / nrow(trilemma_db)) * 100
  ) %>%
  filter(!is.na(age_group))  # Exclude any NA groups

print(age_groups)  





